Arbitrary-Shape Dielectric Particles Interacting in the Linearized Poisson–Boltzmann Framework: An Analytical Treatment

This work considers the interaction of two dielectric particles of arbitrary shape immersed into a solvent containing a dissociated salt and assuming that the linearized Poisson–Boltzmann equation holds. We establish a new general spherical re-expansion result which relies neither on the conventional condition that particle radii are small with respect to the characteristic separating distance between particles nor on any symmetry assumption. This is instrumental in calculating suitable expansion coefficients for the electrostatic potential inside and outside the objects and in constructing small-parameter asymptotic expansions for the potential, the total electrostatic energy, and forces in ascending order of Debye screening. This generalizes a recent result for the case of dielectric spheres to particles of arbitrary shape and builds for the first time a rigorous (exact at the Debye–Hückel level) analytical theory of electrostatic interactions of such particles at arbitrary distances. Numerical tests confirm that the proposed theory may also become especially useful in developing a new class of grid-free, fast, highly scalable solvers.


INTRODUCTION
Electrostatics is pervasive in the physical realm. Especially at the nanoscopic and atomistic scales, it rules plenty of relevant phenomena and plays a crucial role, for instance, in biomolecular systems and biomolecular interactions (e.g., protein−protein association 1 ), colloidal solutions, and atmospheric and plasma physics. During biomolecular recognition and binding, electrostatics provides specificity by helping to tune the delicate balance between the desolvation of the binding interfaces and the descreened charged and polar interactions. 1 At larger distances, the long-range nature of electrostatics allows the interacting partners to spend enough time in close proximity so as to increase the probability of finding the best mutual conformation for binding. 2 The continuum electrostatics description falls into the category of implicit solvent models. 3 It represents the system as a low polarizable medium embedded in a highly polarizable one, which may contain dissolved ionic species described as average spatial densities. The Poisson−Boltzmann equation (PBE) and evolutions thereof 4,5 are continuum electrostatics models. The PBE proved to be very valuable in interpreting many biophysical phenomena. It quite accurately describes the electrostatic interactions occurring between charged particles in solution by combining the electrostatics theory of continuum media with a mean-field approach for the electrostatic potential in the solution. 6 In this context, the PBE applied to molecular dynamics snapshots was used to characterize porosity and solvent-mediated interactions in nucleosome core particles. 7 Thanks to its lower computational demand with respect to an atomistic description, the PBE can also be convenient to study supramolecular structures. For example, the electrostatic potential distribution derived via PBE has been used to improve the accuracy of a multiscale generalized Born model, applied to a 40-nucleosome structure, 8 even though the direct use of a PBE solution via traditional grid-based solvers (such as DelPhi 6,9 or APBS 10 ) becomes impractical for very large systems.
Debye and Huckel (DH) proposed 11 a continuum method for the estimation of the solvation free energy of spherical ion in 1923. In this model, the energy arises from the electrostatic interaction between the ion and the mean potential generated by the surrounding counterion cloud and by the polarizable solvent medium. In their approach, the Poisson equation is solved for the electrostatic potential Φ inside the spherical ion while in solution the free charge linearly responds to the local potential. This formalism corresponds to the linearization of the PBE in spherical symmetry 11−14 (let us recall that in general the PBE is a second-order nonlinear elliptic partial differential equation 15 ). The linearized PBE (LPBE) (for the above reasons also called the DH equation) is often used for low charged systems 16 originating sufficiently small potentials Φ (i.e., as < which is equivalent to |Φ| < k B T/e = 25.7 mV at a room temperature of 25°C, where k B , e, and T are the Boltzmann's constant, elementary charge, and absolute temperature, respectively 12,17−21 ). This, however, does not diminish the importance of studying the behavior of the LPBE also in the case of highly charged objects: their electrostatics may still be correctly described at sufficiently long distances (as compared to the Debye length) by the usual DH approximation provided that the sources of the electric field are properly renormalized. 22−32 (See also the recent ref 33 for additional comments concerning the ranges of applicability of the DH theory.) This once again emphasizes the importance of a thorough study of the DH approximations, both theoretically and numerically, and justifies the constant stream of works related to the LPBE (see recent refs 2, 5, and 34−44, and references therein).
Since the initial developments, considerable efforts have been invested into extending DH theory and deriving analytical approximations for the solution of the LPBE in the more complex case of two interacting particles of spherical shape. Recently, Filippov and Derbenev and co-workers have published several works 19,20,45 exploring the electrostatic force between two charged polarizable spheres immersed in an electrolytic solution or in equilibrium plasma. Other relevant studies important to mention here are refs 2, 12, 40, 41, and 46−52. Let us note that the existing literature largely focuses on treating the case of a system of two spherical particles with special symmetries also in the charge (especially, azimuthal symmetry). Notable exceptions to this are ref 47 (aimed at finding simple analytical expressions for the interaction energy by fitting them to a numerical solution of the DH equation in the case of two equal-radii spheres), ref 48, which extends the results of the previous work to two spheres with unequal radii, ref 50 (aimed at treating interactions with the imposed nonuniform surface potentials), ref 51 (estimating the interaction energy for two rigid globular proteins with arbitrary charge distributions at large separations), and ref 53 (electrostatic treatment of two permeable spherical shells). Let us also note the very recent ref 40 that is the first to provide, in November 2020, the general rigorous (exact at the DH level) treatment of many-sphere systems. On the basis of that work, a further assessment of pair-and multisphere effects is made in ref 41. More detailed review of the literature can be found in ref 2. Albeit the problem of analytically describing the interactions between two conducting spheres is quite long-standing, similar studies on the interactions of two polarizable dielectric spheres (for instance, within the DH framework) arose only relatively recently and are still growing. 19,54−60 It is now known that polarization can strongly influence the electrostatic interactions between dielectric particles, especially at close interparticle separations, and lead to rather counterintuitive effects [that go beyond the scope of the standard Coulombic/singly screened (DLVODerjaguin−Landau−Verwey−Overbeek) interaction terms neglecting polarization], such as the attraction between like-charged particles. 2,19,51,54,55,57,61−64 These effects were observed mostly numerically, but also experimentally, and especially for dielectric spheres. Some partial and approximate analytical results toward the quantification of higher-order terms, which go beyond the conventionally used Coulombic/ singly screened ones, for the potential and interaction energy in a two-sphere system were obtained in refs 12, 46, 51 (doubly screened terms), and ref 40 (triply screened terms), while no general results were known for higher screening orders until the very recent ref 2 to the best of our knowledge. Thus, their study and analytical quantification still remain of great importance. In this respect, in the recent ref 2 the authors presented novel two-center spherical re-expansions that are free of any restrictive symmetry assumptions and improve on the previous developments bypassing the conventional expansions in modified Bessel functions. On the basis of them, they constructed asymptotic expansions in ascending order of Debye screening terms for the electrostatic potential and the total electrostatic energy in the case of two spheres bearing arbitrary charge distributions. This made it possible to explicitly quantify all (k ≥ 0)-screened terms of the potential coefficients and electrostatic energy and thereby to refine a number of partial approximate results previously reported in the literature (see details in ref 2, section II C) for any twosphere system. In the same work, it was further demonstrated that even in the (simplest) case of two centrally located point charges the (k ≥ 2)-screened terms may significantly exceed the conventional singly screened (k = 1) DLVO term. 2,19,41,65,66 This imbalance can only increase when higher-order multipoles are present or when particles have a large dissimilarity (in size, charge, etc.). This emphasizes the importance of developing a rigorous (exact at least at the DH level) electrostatic analytical theory for arbitrary-shape polarizable dielectric particles at arbitrary interparticle separations and without any assumption on charge or system symmetries and particle sizes (as often found also in the recent literature see the detailed overview in ref 2, section III). It is worth noting that, apart from the extensively discussed case of spherical particles, recently exact (analytical) results were also obtained for interaction between a charged dielectric sphere and a planar surface 67 and between two dielectric spheroids 21 in the Poisson limit, i.e., at zero ionic strength. Moreover, recent studies concerned the interaction between cylinders and a flat surface 68,69 and the numerical evaluation of forces in a cylinder−sphere system 70 based on the so-called surface element integration approach (proposed in refs 71 and 72), which attempts to extend the DLVO theory to the interaction of differently shaped particles with a (relatively) flat surface. Here, let us note that the DLVO theory was originally developed for treating spherical particles (colloids), while here it approximates a shape by an "equivalent" spherical diameter. 68, 69 Finally, let us also quote a very recent ref 73, which derives a simple closed-form formula for the apparent surface charge and the electric field generated by a molecular charge distribution in aqueous solution (in the Poisson limit, i.e., at zero ionic strength). However, no rigorous general analytical solution for the case of two polarizable arbitraryshape dielectric particles is currently known, to the best of our knowledge. By rigorous and general throughout this paper we mean both exact at the DH level and free from any restricting hypothesis on geometry or symmetry of the system. The current paper aims at bridging this gap, namely: (1) In order to rigorously treat the mutual polarization of arbitrary-shape particles at arbitrary distances we derive a novel spherical re-expansion result for the LPBE solution. In this approach, no restrictive assumptions on either the symmetry of potentials/charge distributions or on the ratio of r i and R (see section 2 below for definitions of all symbols) are imposed. This is an interesting advance with respect to the existing literaturesee, for instance, refs 2, 12, 19−21, 40, 41, 49−53, 74−81, which require r i < R (see details in section 3.1 and Appendix A below).
(2) On this basis, we derive relations 11 and 12 for determining the potential expansion coefficients both inside and outside two arbitrary-shape dielectric particlessee section 4.1 for details. These relations do not rely on any restrictive assumption and lead to known expressions, such as those in the recent ref 21, for the particular case of azimuthally symmetric interactions of two dielectric spheroids at zero ionic strength (the proof of this fact, however, requires some rather fine mathematical calculations which are postponed to Appendix F).
(3) These relations allow us to construct small-parameter (∝ e −κR /R, see details in section 4.2) asymptotic expansions for potential coefficients and the corresponding total electrostatic energy in ascending order of Debye screening, hereby generalizing the results of recent ref 2 (see section 4.3).
(4) Finally, we perform a brief numerical benchmarking of our analytical theory against the finite-differences based DelPhi ver. V numerical solver 6,9 on several model numerical examples (section 5). Unlike conventional grid-based approaches, our methodology requires no external box boundary conditions and computation time is relatively independent of the distance between particles. Importantly, being grid-free, it does not suffer from numerical artifacts associated with the discretization of the equation. Numerical tests show that the calculation time using the theory proposed in this article can be several orders of magnitude smaller than the corresponding calculation times in DelPhi. Interestingly, different contributions to the potential can be calculated separately with ease. This paper is organized as follows. Section 2 formulates the problem of two interacting dielectric particles relying on the LPBE (DH) model. The transmission conditions treatment and the derivation of novel two-center re-expansion are presented in section 3. Section 4 presents the derivation of relations for determining the potential coefficients and smallparameter expansions. Section 5 demonstrates several numerical tests. Finally, technically subtler derivations, proofs, and auxiliary topics, that are instrumental in (and integral for) this study, are postponed to Appendixes.

ELECTROSTATIC PROBLEM FORMULATION
Let us consider a general system consisting of two nonintersecting dielectric particles i and j, with dielectric constants ε i and ε j . We adopt two spherical coordinate systems with their origins associated with centers x i and x j of the particles (let us note that since the LPBE is a Helmholtz-type equation, it cannot be solved in the standard bispherical coordinate system through separation of variables, see ref 19). Without loss of generality, one can assume that x i and x j lie on the Cartesian axis Z, while the axes X and Y are fixed. The corresponding particle surfaces are then parametrized in these spherical coordinate systems by the radial distances a i and a j depending on the (polar and azimuthal) angles, see Figure 1. The particles are separated by a distance R between their centers. Without loss of generality we will assume henceforth that = i 1, 2 and j = 3 − i. These particles are in an electrolytic solution (for instance, water, and mobile ions) with dielectric constant ε sol and Debye length κ −1 .
The electrostatic potential Φ in,i inside the ith where r i is the radial coordinate of ∈ x  3 measured from the center x i of the ith particle (so that r i = ∥r i ∥, r i = x − x i ), and ρ i (x) denotes the charge density inside the ith particle. Simultaneously, consistent with the Debye−Huckel (DH) model, the potential Φ out,i in the surrounding medium caused by the presence of the ith particle satisfies the LPBE: 12,[19][20][21]46 Due to the superposition principle the self-consistent total electrostatic potential Φ(x) of the whole system is then 12,[19][20][21]46 Ä = i 1, 2, subjected to the following boundary conditions on the particles' surfaces: where n i is the unit normal vector and σ i is a permanent free charge density distribution on the surface r i = a i of the ith particle (if any); since we are further interested in dielectric systems with no fixed free surface charge we can assume σ i = 0 (with no loss of generality of considerations, since formulations with and without fixed surface charge are essentially equivalent from the mathematical point of view, at least for the spheres; see refs 40 and 83). The notation r i → a i ± here indicates approaching the surface of the particle from the interior (−) or the exterior (+) side. Also, A ≔ B or B ≕ A denotes that the value of A is determined (defined) by the value of B.
The general solution of eq 1 can be represented in the form where Φ̂i n,i is the given particular solution to 1 that represents the standard Coulombic potential in infinite space for the distribution ρ i (x); in particular, explicit singling out of the Φ̂i n,i term provides a convenient way to extract the self-energy contributions from the total electrostatic energy of a system (see, e.g., refs 2, 12, 40, 82). Then, introducing dimensionless radial coordinates rĩ ≔ κr i and denoting aĩ ≔ κa i , = i 1, 2, eqs 1 and 2 boil down to the following homogeneous equations: = i 1, 2 and Δ rĩ denotes the Laplace operator with rĩ as the radial spherical coordinate.
Physically feasible general solutions to 6 (such that |Φ̃i n,i | < ∞ as r̃i → 0 + and Φ out,i → 0 as r̃i → + ∞) can be expanded in modified Bessel functions of the second kind, K n+1/2 (r̃i) (Macdonald functions) 84 and associated Legendre polynomials is the nth standard Legendre polynomial) 85 in the real-valued form as follows: with some real-valued expansion coefficients L nm,i , M nm,i , G nm,i , H nm,i to be determined from boundary conditions 4. Appendix I briefly summarizes the minimal necessary information on the modified Bessel functions used in the text. Let us also note, however, that many authors 40,41,51,78 prefer to express potentials 7 in terms of complex-valued spherical harmonics instead of using the real-valued ones (that is, cos(mφ) P n m (μ i ), sin(mφ) P n m (μ i )); this case and the corresponding re-expansion 65 for the DH potential are discussed in Appendix G.
Finally, let us briefly recall that the total electrostatic energy (within the LPBE framework) is given by 6 where ρ fixed is the fixed charge distribution (of any kind, see 1) present in the system. Energy = R ( ) of a given two-particle configuration ( Figure 1) can also be decomposed as 2  The main difficulty in determining expansion coefficients in 7a and 7b from the boundary conditions 4 is that the expansions for Φ out,i (rĩ , θ i , φ) and Φ out,j (rj , θ j , φ) refer to different spherical coordinate systems and corresponding spherical harmonics. For instance, in order to impose boundary conditions 4, the authors of recent refs 19−21 propose to reexpand the potential, say Φ out,j , in terms of coordinates (and corresponding orthogonal Legendre polynomials) of the other sphere i; let us note that this is quite a conventional approach which is followed by many authors, see refs 2, 12, 19−21, 40, 41, 49−53, 74−81, allowing one to handle the corresponding boundary conditions correctly from the mathematical point of view. Let us also note that, in contrast to the well-known works in refs 12, 49, and 74, the theory built in refs 19−21 does not make use of the additional reflection symmetry about the plane bisecting the line connecting the spheres' centers and the corresponding equality of the expansion coefficients of Φ out,i and Φ out,j , which rely on the assumption that the radii of the spheres are equal. Thus, the expansion of Φ out built in refs 19−21 is in principle applicable to the case of spheres with different radii. However, the theory and re-expansions presented in refs 19−21 assume the azimuthal symmetry for the potentials 7a and 7b (i.e., independence of φ) and therefore would not be able, e.g., to deal with an arbitrary orientation of the free dipoles located inside the dielectric spheres. Here, we intend to fill this gap of refs 19−21 and to expand upon the corresponding re-expansions to include general cases devoid of any angular symmetry. Another important feature of the re-expansion presented here is that it does not impose the restrictive inequality r i < R, in contrast to re-expansions derived in the existing literature. 2,12,19−21,40,41,49−53,74−81 This allows us to consider the case of very close arrangement of arbitrary highly irregular dielectric particles. As an example, Figure 2 illustrates the simplest example of such a situation when two flat thick circular dielectric disks are located very close to each other (let us note that despite the considerable interest to the interactions between two flat membranes/disks 14,18,86,87 we are not aware of any complete, DH-exact, analytical description of this kind of system).
Namely, to this end, we advance the following representation (re-expansion) of the potential Φ out,j : where the re-expansion coefficients b nml are determined via 29 and R̃≔ κR, = i 1, 2, j = 3 − i. The quite technical derivation of 9 is given in Appendix A.
Let us emphasize that 9 provides an expansion of the potential Φ out,j (originally referred to coordinates of the jth spherical system and having harmonic expansion coefficients G lm,j , H lm,j ) through harmonics and coordinates referenced now to the ith system. Analytical properties (alongside some important particular cases) of the re-expansion coefficients b nml are described in Appendices A and B.
3.2. Numerical Calculation of Potentials in Practice: Truncating the Re-expansions and Approximating the Re-expansion Coefficients. Since the quantities (∑ l = m +∞ b nml (rĩ , R̃)G lm,j ) and (∑ l = m +∞ b nml (rĩ , R̃)H lm,j ) in 9 as well as the re-expansion coefficients b nml (r, R̃) of 29 in general contain infinite sums, in practical calculations, to determine the potential expansion coefficients L nm,i , M nm,i , G nm,i , H nm,i , = i 1, 2, one needs to apply a truncation to a finite number of terms. This is usually done according to the required accuracy, which is often estimated by tracking the evolution of some key quantity, e.g., the electrostatic energy of the system as done in the convergence estimate for grid-based solvers. 6,9 Interestingly, only the energy components depending on potential coefficients subjected to further changes need to be recalculated; see 7a and 8. To this end, we propose and then numerically benchmark (section 5) the approximation methodology, the simpler "azimuthally symmetric" version of which for m = 0 was proposed in ref 19 and successfully verified to be effective in refs 20 and 21. The approximation methodology we propose consists of the following two points: (1) only coefficients b nml (r, R̃) with n + l − m ≤ n max are to be calculated, while all of the others are assumed to be zero; (2) further additional constraints s ≤ n max − l and k ≤ n max + m − l − s are enforced on the infinite series 29, where n max ≥ m, is a given fixed user-defined threshold. For m = 0, the proposed approximation methodology simply boils down to that of refs 19 and 20. Simple algebraic calculations indicate that this approximation methodology calculates the exact values of the coefficients b nmm (r, R̃) (see 36) if n max ≥ n is taken; however, this is not the case, e.g., for the coefficient b nml (r, R̃) with general triplet (n, m, l) of indices (an illustrative example of the convergence of the re-expansion coefficients, approximated by the methodology just described, is given in Appendix B.3).

EXPANSION COEFFICIENTS FOR THE POTENTIALS,
SMALL-PARAMETER EXPANSIONS: THEORY 4.1. Derivation of Relations Governing the Potential Coefficients. Determining the unknown potential expansion coefficients L nm,i , M nm,i , G nm,i , H nm,i of 7 completely solves the problem of finding the electrostatic potential.
With using 5, 7, and 9, after algebraic transformations boundary condition 4a acquires the following expanded form: i n , i i where, to shorten the recording of formulas, it is denoted that so that the entire surface of particle i is covered), = i 1, 2, j = 3 − i, and ∑ n m , denotes the double sum (likewise to 7) over indices (n, m) with n ≥ m ≥ 0 or n ≥ m ≥ 1 for the expressions involving coefficients L nm , G nm , or M nm , H nm , respectively. Then, multiplying both sides of this equality by  (10) and integrating over Ω i ≔ {(θ i , φ)|0 ≤ θ i ≤ π, 0 ≤ φ < 2π} (i.e., over the entire surface of particle i) with weight sin θ i , one gets the following linear systems with respect to the unknown coefficients of 7: ;cos   note that in the (simplest) case of a spherical surface (i.e., if a i (θ i , φ) = constant independent of angles θ i , φ), functions 10 constitute a complete orthogonal set on a sphere parametrized by Ω i so that the integral ∬ Ω i (·) sin θ i dθ i dφ of the product of arbitrary two such functions with indices (n′,m′) and (n″,m″) is zero if (n′,m′) ≠ (n″,m″) (in particular, complex-valued spherical harmonics Y nm (θ,φ) are constructed from this basis and fulfill the same orthogonality relation 82 ); see 50 below. Then, systems 11 boil down to the identities just resulting from simple collecting/equating the expansion coefficients (of all the functions involved in 4a) at Fourier spherical harmonics sin(mφ)P n m (μ i ), cos(mφ)P n m (μ i ) of the same orderssee, e.g., relation 22 below. This special case is discussed in more detail later in section 4.3. In the general case, when a i (θ i , φ) does not describe a sphere, the above-described approach of treating boundary conditions essentially follows the idea of spectral Galerkin residual orthogonalization procedure, 88−91 with the trial and test functions spaces being spanned by set 10.
Let us account for the second boundary condition, that is, relation 4b. Following the same approach as in the previous case of boundary condition 4a and using expressions 46 and 47 to treat differential operators n i · ∇, and relations k x ( ) , eq 8.731), we arrive at the following linear systems with respect to the unknown coefficients of 7: ;cos where coefficients g n′m′,nm i; cos , h n′m′,nm i; cos , i n′m′,nm i; cos , j n′m′,nm i; cos , k n′m′,nm i; cos , l n′m′,nm i; cos , and n n′m′ i; cos are given by 49. The values of g n′m′,nm i; sin , h n′m′,nm i; sin , i n′m′,nm i; sin , j n′m′,nm i; sin , k n′m′,nm i; sin , l n′m′,nm i; sin and n n′m′ i; sin are defined in the same way, except that the integrals 49 use sin(m′φ) instead of cos(m′φ) in their integrands.
General systems 11 and 12 read in matrix form as follows: ; ; (13a) ; ; (13b) (13c) i; cos }, and corresponding matrices/vectors with subscript i; s are defined by coefficients with superscript i; sin in exactly the same way. Let us also note that, by construction, indices (n′,m′) enumerate the rows in (sub)matrices 13, while (n,m) enumerate their columns so that (n′,m′) run 0 ≤ m′ ≤ n′ in matrices/vectors with subscript i; c and 1 ≤ m′ ≤ n′ in those with subscript i; s, while (n,m) run 0 ≤ m ≤ n in matrices A, C, E, G, I, K, and 1 ≤ m ≤ n in B, D, F, H, J, L; in addition, when using approximation approach of section 3.2, the corresponding indices n′ and n are bounded from above by n max . Next, mutually swapping indices i and j in 13 one also gets the similar four-equation system but with i and j interchanged. Thus, combining 13 and the corresponding system with i and j interchanged, we assemble the following global linear system with the block matrix composed of separate matrix blocks (submatrices) and the (global) unknown column-vector composed of separate column-vectors i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z z (there 0 denotes a matrix with all-zero entries). Solving the global linear system 14 one An alternative way for determining these coefficients (yielding explicit expressions for them as well as also allowing us to construct the corresponding small-parameter expansionssee the next section 4.2) is discussed in Appendix E.
Let us note that in the case of a configuration with only one (say, ith) dielectric particle present, system 14 simplifies to i k j j j j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j y { z z z z z z z z z z z z z z z z z i k j j j j j j j j j j j j j j j j j j j j j j j y = i 1, 2, where addends with superscript (k) are k-screened, that is, tend to zero as (e −κR /R) k as R → + ∞ (or, to be mathematically tidier, they are of order O((e −κR /R) k ) as R → + ∞).
where addends with superscript (k) are k-screened. Note also that expansions for ̂̂i and ̂̂i contain only even-screened and odd-screened addends, respectively; in particular, doing so one can obtain that  where O(e −κR /R)-terms have been underlined (just for better traceability the screening origins of the corresponding lefthand sides). Thence, using Neumann matrix series and employing 17 in 59 one can finally get the desired addends of expansion 16 for G i ranging in ascending order of Debye screening; in particular, by doing this one can write down the first few addends: } of expansion 16 of G i further allow us to construct expansions for the remaining vectors H i , L i , and M i . Indeed, employing Neumann matrix series and identity 56 one has 1 ; , or, by collecting terms of the same order of screening, we conclude that Then, using 54 one finally obtains In turn, employing the coefficients 20 in 7a and further using 8, 5, and 3 immediately yield 2 the k-screened components k ( ) of asymptotic expansion (16) , 0 which coincides with ref 2, eq 43; note that integrals 48 containing sin(mφ) completely nullify in 50, for instance, b n′m′,nm i; cos , d n′m′,nm i; cos and f n′m′,nm i; cos involved in 11a, whereas only those containing cos(mφ) with m = m′ survive. Using 50 in 54 we further arrive at 51, employing which in 55 we obtain 52. Matrix elements 52, in turn, lead to 53 which coincides with ref 2, eq 45, governing the external potential coefficients in the problem of interaction of two spheres. Furthermore, using 50, 51, and 52 after some algebraic transformations it can be shown that asymptotic expansions addends 18, 19 and 20 for the potential coefficients boil down to those of ref 2, eqs 60 and 69, for the case of a two-sphere system (which in turn were also extensively tested numerically for convergence behavior 2 and were shown to generalize a number of approximate results previously reported in the literature, see ref 2, sections V and VI), e.g., in the simplest case of two centrally located point charges q i and q j (thus, Φ̂= immediately follow. 2 In the particular case of equally sized spheres a i = a j = a the latter relation reduces to the well-known DLVO interaction energy 18,63 0 sol 2 ; see ref 2 concerning the further details on the full analytical solution to the problem of two interacting spheres bearing arbitrary multipoles. Thus, the rigorous (exact at the DH level) electrostatic theory and the corresponding asymptotic small-parameter expansions built in section 4.2 indeed naturally generalize those of ref 2 for the two-sphere system to a much more complicated case of two arbitraryshape interacting particles.

NUMERICAL MODELING AND RESULTS
In this section we present the numerical illustration of the developed theory; we also compare the corresponding numerical results with those provided by the DelPhi LPBE solverit is especially well-suited for our current purposes since the version of DelPhi introduced in ref 9 has built-in support for the insertion of the simple geometric objects (e.g., dielectric particles in the form of spheres, cylinders, cones, etc.).
Let us note that the advantage of the approach proposed in the current paper is the relatively small size of the linear system to be solved for calculating the electrostatic potentials in contrast to conventional grid-based methods (where variables/ degrees of freedom of the corresponding linear system are usually associated with grid points, and accurate solution of the problem inevitably entails the usage of fine grids), e.g., the matrix of 14 is N × N with N = 4(n max + 1) 2 and, as we will observe from the numerical experiments of the current section, relatively small n max 's are usually sufficient to produce quite accurate solutions.
Typical values used in the calculations are ε i = ε j = 2, ε = 80 sol , and κ −1 = 8.071 Å (these solvent parameters are quite typical for systems of biophysical interest, representing an aqueous solution with 0.145 M physiological NaCl concentration at a room temperature of 25°C). The DelPhi's settings 6,9 perfil (the percentage of immersion of the physical system into the computational cubic box, to be minimized in order to make the approximation in the external boundary conditions more acceptable) and scale (the reciprocal of one grid spacing, grids/angstrom, to be maximized in order to obtain a finer mesh resolution), which directly affect the time and accuracy of numerical calculations, are described below. Finally, as a stopping criterion in DelPhi we set the value of parameter maxc (grid potential maximum change threshold, in k B T/e units, where k B , e, and T are Boltzmann's constant, elementary charge, and absolute temperature, respectively) to 10 −7 (unless explicitly stated otherwise).
Numerical experiments are discussed in the next section, 5.1. In addition, in section 5.2 we also propose and discuss a simple approach, based on the Tikhonov regularization theory, to improve the accuracy and stability of potential/energy calculations at high n max .
where 2z 0 is the cone height, z 0 = 10 Å, and two point charges q 1 = q 2 = 10e are situated in their centers x 1 , x 2 . The corresponding graphical representations for the distances R = 50 Å and R = 25 Å are shown in Figure 3. These cases are The Journal of Physical Chemistry B pubs.acs.org/JPCB Article azimuthally symmetric (thus, potential coefficients with m > 0 nullify and therefore do not need to be calculated, in general).
Then, Figure 4 demonstrates the total electrostatic energy for R = 25 Å and R = 50 Å at different n max 's calculated by the proposed methodology; the figure shows that stabilizes rather quickly and practically stops changing after n max = 9. Table 1 reports in more detail the numerical values of at R = 50 Å for some n max 's together with the corresponding calculation timings in serial and parallel (values in parentheses) modes. Table 1 clearly demonstrates that the calculation time using the theory proposed in this article can be several orders of magnitude less than the corresponding calculation times in DelPhi (see Figure 5 for more details on the calculation times using the proposed theory in the current example, in both cases of serial and parallel computations). Figure 5 shows that the computation time in our approach has little dependence on R (moreover, since only the matrices E i;c , E i;s , F i;c , F i;s , K i;c , K i;s , L i;c , and L i;s depend on R, and this dependence is of the order of O(e −κR /R), which causes these matrices to decay rapidly, see section 4.2, then in practice, the integrals 48 and 49 may be calculated even faster for larger R within a given accuracy). In contrast to this, grid-based approaches must cover the whole system under interest by a computational domain; thus, decreasing R decreases the grid size too and, respectively, the total time of calculations, e.g., for R = 25 Å the corresponding DelPhi calculation times (for the same scale and perfil as consecutively listed in Table 1)     a Times indicated in parentheses refer to parallel computing on four parallel workers (see more details in Figure 5).  (Figures 6a and 7a) are visually indistinguishable from those calculated by this approach (see Figures 6b and 7b) at n max = 6. The corresponding absolute errors max Ω h |ϕ − ϕ DelPhi | over the whole Ω h for n max = 6 are 0.28 and 0.20 for R = 25 Å and R = 50 Å, respectively (an increase in n max makes it possible to further reduce these errors, e.g., at n max = 10 they become equal to 0.18 and 0.14, respectively). All integrals 48 and 49 were calculated here using the builtin MATLAB functions (integral/quadgk) with the default MATLAB's accuracy settings (AbsTol = 10 −10 , RelTol = 10 −6 ). Let us also note that the step in n max simply implies adding new rows and columns to matrices A i;c , (the similar identity for a 2 can be easily derived by inverting here θ 1 to π − θ 1 ), where 2z 0 is the cylinder height, z 0 = 10 Å, h is the cylinder radius, h = 50 Å, and two point charges q 1 = q 2 = 10e are situated in their centers x 1 , x 2 . As in the previously considered case of a system of two cones (see Figure 4), for sufficiently large n max , energy stabilizes and practically stops changing. For instance, in the case shown in Figure 2a Table  2; it can be seen that already at sufficiently small n max equal to 6 and 10, the proposed approach provides values very close to that of Delphi (relative error about 0.01%).    The experience on grid-based PB solvers taught us to perform energy calculations while preserving the relative position of the system fixed with respect to the grid. As shown in the very simple case of two approaching charged amino acids, see for instance Figure 11 in ref 93, the calculation of the potential of mean force is particularly delicate in this respect since by construction it does not meet the requirement of fixed geometry. Here, we show, as a proof of concept, arginine and glutamate charge distributions placed into two cylindrical dielectric particles (see section 5.1.2 for the corresponding surface parametrizations). We assume that x 1 and x 2 coincide with the corresponding geometric cylinders' centers and x 1 = (0, 0, 0) (glutamate charge distribution) is fixed while x 2 = (0, 0, R) (arginine charge distribution) changes from R = 8.54 Å to R = 14.94 Å with a step of 0.1 Å (see Figure 10). Then, Figure 11a shows the results of the calculation of the total electrostatic energy with the proposed methodology and using DelPhi with the parameters perfil = 80, scale = 3, maxc = 10 −5 , which correspond to those used in ref 93; in addition, the results for scale = 15 (which is, in general, an incredibly large value as compared to scales usually used in DelPhi for calculations in systems of biological interest 6,9,93 ) and maxc = 10 −7 are also shown. One can observe that the spurious oscillations pollute the energy profile calculated by DelPhithese are caused by the numerical (grid) artifacts that are due to the discretization of the equation. In a more realistic case, where the molecular surface would be used rather than a basic geometric model surface, the smallest atomic radius for arginine in the CHARMM parameter set, 93 which is about 0.22 Å for some polar hydrogens, would even more largely contribute to this phenomenon and require extremely, and practically unfeasible, fine grids. 93 At the same time, Figure 11a shows that the energy profile calculated by the approach proposed in the current paper is free from such shortcomings. Finally, Figure  11b illustrates the convergence of as n max increases.    21, which builds the rigorous theory of electrostatic interactions of two spheroids in the azimuthally symmetric case at κ = 0, observed in their calculations that the corresponding linear systems governing the potential coefficients (in our case, system 14) may be illconditioned for large n max possibly leading to numerical instabilities/artifacts and thence to a loosening of numerical calculations with a further increase in errors in the potential as n max increases. For instance, Figure 12 illustrates how the 2norm condition number (CN) (i.e., the ratio of the largest singular value to the smallest one) of the linear system governing the potential coefficients grows with increasing n max in the example with two cones considered above. Let us note that, although CN is a rather rough characteristic, it can still serve as a general indicator of how sensitive a linear system is to numerical errors (inaccuracies that arise when calculating the coefficients of the system, round-off errors in the numerical solution of the system itself, etc.) and how these errors can affect its solution. 94 Large values of CN indicate that the numerical solution results may be inaccurate/unreliable (e.g., when CN ≳ 10 16 and double-precision floating-point arithmetic is used). Unfortunately, no specific solution for illconditioning was provided in ref 21. Thus, to enhance the stability and robustness of calculations we adopt the following simple approach (which conceptually follows the Tikhonov regularization theory 95 ): namely, instead of solving the original (possibly ill-conditioned) system 14 represented here in the matrix form Ax⃗ = b ⃗ we solve the perturbed system (A′A + αE)x⃗ = A′b ⃗ , where A′ denotes the conjugate transpose of A, α > 0 is a regularization parameter (so that x⃗ depends on α now: x⃗ = x⃗ α ), and E is some (wellconditioned) symmetric positive-definite regularizer (e.g., the simplest choice for E, also tested/used in our numerical experiments, is just the identity matrix). Since the perturbed matrix A′A + αE is symmetric, positive-definite, and wellconditioned (thanks to the regularizing addend αE), the corresponding perturbed system can now be effectively handled using, e.g., Cholesky or LDL decompositions. 94 For the choice of α we followed the idea of the so-called noise level-free quasi-optimality criterion 95−97 employing the geo-metric sequence α = α i = α 0 q i (where 0 < q < 1, i = 1, ..., M) and then selecting α i which gives the smallest discrepancy in the 2-norm (or, alternatively, one may also rely on the values of the classical penalized least-squares functional instead; however, in our numerical experiments of section 5 this led to almost the same results). Figure 13 demonstrates the results of such a regularization (we have used α 0 = 0.8 5 , q = 0.8, and M = 100 in our numerical experiments) in the problem of two interacting cones (see section 5.1.1) at R = 50 Åit can be seen that the regularized solution behaves in a more stable way (since it is better conditioned). One can also observe from Figure 13 that increasing the accuracy of calculation of the integrals 48 and 49 forming system 14 (as compared to the default accuracy of built-in MATLAB integration functions) as expected improves the numerical solution. We also note that the described regularization procedure could be obviously applied regardless of n max ; however, as our numerical experiments suggest, regularization begins to play a role only for sufficiently large n max (e.g., in our numerical experiments, for n max > 20; see Figure 13) while for smaller n max the regularized and nonregularized solutions practically coincide. At the same time, the computational cost for such a simple regularization is negligibly small compared to the overall time of calculating the integrals 48 and 49 (especially thanks to the small size of the linear system governing the potential coefficientssee the comments on this at the beginning of section 5).
Let us finally note that despite that the numerical solution converges rather rapidly with increasing n max (as we can observe from the numerical experiments of section 5.1) and such a simple regularization methodology, considered in the current subsection, usually significantly enhances the stability/ reliability of calculations and alleviates the process of numerical solution overall, the (ill-)conditioning of the linear systems governing the potential coefficients may still be a bottleneck  The Journal of Physical Chemistry B pubs.acs.org/JPCB Article issue in the practical computational/numerical applications of the proposed approach and thus needs to be properly addressed in the future studies (beyond the current proof-ofconcept analytical work). In this respect, we foresee at least the following two possibilities: (1) investigating better choices for the regularizer E and regularizing parameters {α i } in the current regularization scheme as well as adopting other techniques and rules (see refs 95−98) for estimating and regularizing the potential solutions; (2) developing ad hoc (i.e., specialized for system 14) preconditioners with subsequent usage of iterative methods for solving 14. However, we do not have the possibility to pursue these directions further here.

DISCUSSION AND CONCLUSIONS
This paper considers the interaction of two arbitrary-shape polarizable dielectric particles immersed into solvent assuming that the linearized Poisson−Boltzmann equation holds. In order to rigorously treat the mutual polarization of arbitrary-shape particles at arbitrary distances R, in section 3 we present a novel spherical re-expansion for the LPBE solution. Advancing what can be found in the existing literature (refs 2, 12, 19−21, 40, 41, 49−53, 74−81, 99, 100) neither assumptions on the symmetry of potentials or charge distributions nor on the ratio of r i /R are made. Although the obtained general re-expansion coefficients 29 contain infinite series (in contrast to the r i < R case, where these expressions boil down to more compact finite sums 31 and 32), in section 3.2 we propose and discuss an efficient approximation procedure and validate it numerically (Appendix B.3 and section 5).
On this basis, in section 4 we then derive relations governing the potential coefficients (section 4.1). In turn, they then allow us to construct small-parameter (∝ (e −κR /R) k ) asymptotic expansions for the potential coefficients and for the total electrostatic energy in ascending order of Debye screening (section 4.2). These generalize the results established in recent ref 2 for the case of two spherical particles.
Finally, in section 5, we perform the numerical benchmarking of this analytical derivation validating it against the wellknown grid-based DelPhi numerical solver 6,9 on several model examples. Computational examples have been provided with basic shapes, such as cones and cylinders, which can approximate more complex structures at the nanoscale (the general theory built in the article is suitable for arbitrary-shape particles having an analytical representation a i (θ i , φ), see section 2). Advantages of this approach with respect to conventional grid-based techniques reside in the fact that (i) it is inherently consistent with null boundary conditions for the potential at infinite distance from the solute(s), (ii) its performance is practically independent of the distance R between the particles, (iii) being grid-free it is not subjected to numerical artifacts associated with the LPBE discretization or to the presence of the so-called self-energy, 2,6,82 and (iv) finally, if needed, specific contributions, such as the components arising specifically from the polarization charge at the particle boundaries or ionic contributions can be singled out and studied analytically. 2 Numerical tests show that the calculation time using the theory proposed in this article can be several orders of magnitude smaller than the corresponding ones in DelPhi. Moreover, a simple parallelization scheme, acting on the assembling process of the elements of system 14, which are governing the potentials, can also be appliedsee Figure 5 and corresponding explanations in section 5.
Applications of this theory range from a better way of benchmarking numerical grid based approaches for the LPBE, as well as for a better approximation of their boundary conditions, 101 to allowing a careful study of how geometry impacts on interaction energy, to the treatment of mesoscale systems, approximated as simpler spheroidal or ellipsoidal particles, and their mixtures, 102 in the fields of both biomolecular modeling, supramolecular assemblies, and colloids. Due to the absence of grid artifacts, this approach appears particularly useful for applications such as the calculation of the potential of mean force, where the same relative position between each of the two particles and the grid can hardly be preserved while the relative distance is changed. Current work is ongoing to instantiate the present formalism in the case of conventional atomistic description of biomolecular systems, such as implementing the various definitions of the protein molecular surface. 93 ■ A. DERIVATION THE PROPOSED NOVEL

RE-EXPANSIONS AND THE CORRESPONDING RE-EXPANSION COEFFICIENTS FOR DH POTENTIALS
A.1 Derivation of Re-expansion 9 Taking into account that where R̃≔ κR, = i 1, 2, j = 3 − i, we will use the following Macdonald (note that (−1)!! = 1 is always assumed 92 ). Further, in order to decompose the product rj l P l m (μ j ) we also rely on the following decomposition theorem for the associated Let us comment that the original identity is established in ref 85 for the angle π−θ i in the argument of P k m on the right-hand side, but by using P k m (−x) = (−1) k+m P k m (x) one can easily obtain 26.
The Journal of Physical Chemistry B pubs.acs.org/JPCB Article Finally, we take advantage of the following representation for the product of a Gegenbauer polynomial with an associated Legendre polynomial: where the numerical coefficients h kmls,n exist and are uniquely determined for arbitrary given polynomials P k m (x), C s l+1/2 (x). The representation given in 27, although quite simple, seems to be unnoticed earlier in the literature, and therefore its proof is given in Appendix A.2 alongside with explicit closed-form relations for calculating h kmls,n . Now using 25, 26, and 27, let us now transform the solution Φ out,j (rj , θ j , φ) given by 7b: (let us note that h kmls,n = 0 if s + k − n is oddsee Appendix H) and performing transformations of the sums, we finally arrive at representation 9. The double sum in 29 may also be left unbounded if it is understood that all h kmls,n with n > s + k as well as when s + k − n being odd vanish: Note that, by definition, inequalities n ≥ m and l ≥ m always hold for the indices of b nml (r, R̃).
Let us also note that re-expansion coefficients 29 generalize the corresponding "azimuthally symmetric" coefficients of refs 19 and 20 (see ref 19, eq 11) in the sense that they recover the latter at m = 0 (i.e., in the case of azimuthal symmetry) and r i < R; at the same time, the re-expansion of the azimuthally symmetric potential Φ out,j (rj , μ j ) previously derived in refs 19 and 20 is the particular case of the representation 9 (upon the assumption of the azimuthal symmetry and r i < R conditions).
Re-expansion coefficients 29 also generalize (recover at m = 0) the "azimuthally symmetric" coefficients of refs 99 and 100. Besides, in the particular case of r i < R, coefficients 29 boil down to equalities previously derived in ref 2, eqs 37 and 40: The advantage of relations 31 and 32 is that they do not contain any infinite series of modified Bessel functions in their R-dependent parts; these relations were derived in ref 2 in a completely different way as compared to the current paper (in ref 2 we have employed Gegenbauer−Sonine identities and the generalized Neumann transforms; however, such a method of proof essentially relied on the condition of r i < R, in contrast to what is done in the current paper). In Appendix B.1 we provide an analytical derivation of exact explicit (closed-form) expressions for some typical values of re-expansion coefficients and compare them with 31, which is also used in testing the proposed numerical methodology for calculating b nml (r, R̃) (see section 3.2). Finally, let us also note that, in general, the definition 29, in contrast to 31 with 32, does not represent the expansion of b nml (r, R̃) in powers of R (due to the analytical representation 68a for K n+1/2 (R̃)). Nevertheless, using 70 in 28 immediately yields as r̃→ 0 and R̃→ 0, r̃< R̃, n ≥ 0; thus, for small values of r̃and R̃(e.g., in the weak screening regime as κ → 0) one has the relation Ω̃+ r R ( , ) n 1/2 ≈ π̃+ r R 2 n n 1 . Employing this approximation for Ω l+s+1/2 (r, R̃), the proposed methodology (see section 3.2) allows us to list in 29 all the terms up to (and including) R −(n max +m+1) .

A.2 Explicit Construction of the Coefficients h kmls,n
We start with the known unique representation of an arbitrary monomial power x j through Legendre polynomials P n (x) (see Next, since C s l+1/2 (x) is a polynomial of degree s and P k m (x) = is a polynomial of degree k − m, then their product Thus, the explicit closed-form construction provided in 35 completes the proof of representation 27. Also, additional identities for calculating the coefficients h kmls,n are provided in Appendix H. where δ is a Kronecker delta. Employing this in 29, we immediately arrive at the following:  (36) where s̃= min(r, R̃), S̃= max(r, R̃). We also take advantage of the following exact equality valid for arbitrary n ≥ m ≥ 0: (2) Determination of h m+1,m,m+1,s,n : by virtue of 27 we have P m+1 m (x)C s m+1 + 1/2 (x) = ∑ n = m m + s + 1 h m+1,m,m+1,s,n P n m (x), or using These series follow directly from 25 (when r̃≔ rĩ , l = 0 and μ i = cos θ i = ± 1 there) and are absolutely convergent. Now using 45 in 44 and performing algebraic transformations we arrive at the desired identity 37.

B.3 Numerical Validation of the Re-expansion Coefficients Convergence
As an example of the methodology of section 3.2 for approximation of the re-expansion coefficients, Figure 14 provides the corresponding illustrations of the convergence of the calculated values of b 001 (r, R̃) and b 101 (r, R̃) to the exact values given by 38 and 39. Let us also comment on Figure 14 that, owing to the fact that h kmls,n = 0 if s + n + k is odd (see  , result in the same b 101 (r, R̃). One can observe in Figure 14 the fast convergence of the re-expansion coefficients approximated by the described methodology to the exact ones. One can also observe in this figure that the discrepancy is larger as R̃gets smaller; a further study of the finer convergence properties of the re-expansion coefficients in the vicinity of R̃= r̃remains to be addressed in future work.
The normalized (unit) outward normal vector n i to the surface r i = a i (θ i , φ) of the ith particle can be represented as If Φ̂i n,i (rĩ , θ i , φ) is representable through multipoles 21, the corresponding m n′m′ i; cos takes the following value:    (49) If Φ̂i n,i (rĩ , θ i , φ) is representable through multipoles 21, the corresponding n n′m′ i; cos takes the following value: where we have denoted auxiliary matrices   and vectors Relation 57 and that of with indices i and j interchanged immediately yield Now identities 59, 56, and 54 provide explicit (merely resolved in terms of the coefficients of system 14) expressions for Let us finally note that, alternatively, instead of obtaining 56, 57, and 58, one may also express G i via H i and H j from 55a, and then substitute it into 55b to obtain an equation relating solely H i and H j (completely similar to 57); however, we do not pursue this direction further here.
Indeed, the right-hand side of 60, with numeric coefficients G nm,i (lapl) and H nm,i (lapl) arbitrarily given, represents a general formal solution to the Laplace equation in the ith spherical coordinates (r i , θ i , φ) vanishing at infinity. 84,85 The where r i < R, r j = μ + − R r Rr 2 i i i 2 2 , μ j = (R − r i μ i )/r j , i ∈ {1,2}, j = 3 − i, and with the corresponding infinite series being absolutely and uniformly convergent. Employing 61 in 60 one easily gets the desired re-expansion: Let us emphasize that 62 is obtained independently of the PBE re-expansion theory built in Appendix A and provides the reexpansion of the Laplace solution of general form 60 with arbitrarily given numeric coefficients {G nm,i (lapl) } and {H nm,i }. As well, multiplying expansion 61 by P n m , denoting z = r i /R < 1 there, integrating term by term (it is permitted because 61 is uniformly convergent) and using the orthogonality property which immediately follows from the orthogonality property 66 recast to 27 might be easier for practical calculations, since for given particular values of the indices the last integral is readily computable analytically by any computer algebra system. However, to the best of our knowledge, the integral 67 is unavailable in the literature, so let us find the explicit expression for 67 in a closed form. Due to the oddness of the integrand in 67, h kmls,n = 0 if s + n + k is odd. Let us consider the situation when s + n + k is even. where for the summation |n − k| ≤ p ≤ n + k, p ≥ 2m, p + n + k is even, and C m,m,2m n,k,p , C 0,0,0 n,k,p are Clebsch−Gordan coefficients. 111 It is worth noting that using the Wigner 3-j symbols 111   and 4 F 3 is the generalized hypergeometric function 92,112 of unit argument that is available for calculation by the built-in procedures in most computer algebra systems. Moreover, it was shown in ref 112 that the corresponding hypergeometric series is a terminating one if l, m, k, n are non-negative integers and Re((m + n)/2 − r) > 0, which is valid in our case. ■ I. MODIFIED BESSEL FUNCTIONS K n+1/2 (·) and I n+1/2 (·) Functions K n+1/2 (·) and I n+1/2 (·) of semi-integer order n + 1/2 have the following exact analytic representations: 84 ( 1) e ( )